// Sparse neural network using multiple width 4 neural layer with
// forward connected weight switching as the non-linearity. Uses 
// the fast Walsh Hadamard transform (WHT) as a connectionist device 
// (synthetic fully connected layer.)
// Uses a sub-random pattern of sign flips to prevent the first
// WHT producing a biased spectral response to natural data.
// Uses the P5 (Processing) JavaScript library.
// Sean O'Connor 18 September 2021 
class SWNet4 {
  constructor(vecLen, depth) {
    this.vecLen = vecLen;
    this.depth = depth;
    this.scale = 1.0 / sqrt(vecLen);
    this.params = new Float32Array(8 * vecLen * depth);
    this.flips = new Float32Array(vecLen);
    for (let i = 0; i < vecLen; i++) {
      this.flips[i] = sin(i * 1.895) < 0 ? -this.scale : this.scale;
    }
    for (let i = 0,j=0; i < this.params.length; i += 8) {
      this.params[i+j] = this.scale;
      this.params[i+j + 4] = this.scale;
      j=(j+1)&3;
    }  
  }
  recall(result, input) {
    for (let i = 0; i < this.vecLen; i++) {
      result[i] = input[i] * this.flips[i];
    }
    let paIdx = 0; // parameter index
    for (let i = 0; i < this.depth; i++) {
      wht(result);
      for (let j = 0; j < this.vecLen; j += 4) {
        let x = result[j];
        let a,b,c,d;
        if (x < 0) {
          a = x * this.params[paIdx];
          b = x * this.params[paIdx + 1];
          c = x * this.params[paIdx + 2];
          d = x * this.params[paIdx + 3];
        } else {
          a = x * this.params[paIdx + 4];
          b = x * this.params[paIdx + 5];
          c = x * this.params[paIdx + 6];
          d = x * this.params[paIdx + 7];
        }
        paIdx += 8;
        x = result[j + 1];
        if (x < 0) {
          a += x * this.params[paIdx];
          b += x * this.params[paIdx + 1];
          c += x * this.params[paIdx + 2];
          d += x * this.params[paIdx + 3];
        } else {
          a += x * this.params[paIdx + 4];
          b += x * this.params[paIdx + 5];
          c += x * this.params[paIdx + 6];
          d += x * this.params[paIdx + 7];
        }
        paIdx += 8;
        x = result[j + 2];
        if (x < 0) {
          a += x * this.params[paIdx];
          b += x * this.params[paIdx + 1];
          c += x * this.params[paIdx + 2];
          d += x * this.params[paIdx + 3];
        } else {
          a += x * this.params[paIdx + 4];
          b += x * this.params[paIdx + 5];
          c += x * this.params[paIdx + 6];
          d += x * this.params[paIdx + 7];
        }
        paIdx += 8;
        x = result[j + 3];
        if (x < 0) {
          a += x * this.params[paIdx];
          b += x * this.params[paIdx + 1];
          c += x * this.params[paIdx + 2];
          d += x * this.params[paIdx + 3];
        } else {
          a += x * this.params[paIdx + 4];
          b += x * this.params[paIdx + 5];
          c += x * this.params[paIdx + 6];
          d += x * this.params[paIdx + 7];
        }
        paIdx += 8;
        result[j] = a;
        result[j + 1] = b;
        result[j + 2] = c;
        result[j + 3] = d;
      }
    }
    wht(result);
  }
}

// Fast Walsh Hadamard Transform provide your own scaling
function wht(vec) {
  let n = vec.length;
  let hs = 1;
  while (hs < n) {
    let i = 0;
    while (i < n) {
      const j = i + hs;
      while (i < j) {
        var a = vec[i];
        var b = vec[i + hs];
        vec[i] = a + b;
        vec[i + hs] = a - b;
        i += 1;
      }
      i += hs;
    }
    hs += hs;
  }
}

// Sum of squared difference cost
function costL2(vec, tar) {
  var cost = 0;
  for (var i = 0; i < vec.length; i++) {
    var e = vec[i] - tar[i];
    cost += e * e;
  }
  return cost;
}

class Mutator {
  constructor(size, precision, limit) {
    this.previous = new Float32Array(size);
    this.pIdx = new Int32Array(size);
    this.precision = precision;
    this.limit = limit;
  }
  mutate(vec) {
    for (let i = 0; i < this.previous.length; i++) {
      let rpos = int(random(vec.length));
      let v = vec[rpos];
      this.pIdx[i] = rpos;
      this.previous[i] = v;
      let m = 2 * this.limit * exp(random(-this.precision, 0));
      if (random() < 0.5) m = -m;
      let vm = v + m;
      if (vm >= this.limit) vm = v;
      if (vm <= -this.limit) vm = v;
      vec[rpos] = vm;
    }
  }
  undo(vec) {
    for (let i = this.previous.length - 1; i >= 0; i--) {
      vec[this.pIdx[i]] = this.previous[i];
    }
  }
}

// Test with Lissajous curves
let c1;
let ex = [];
let work = new Float32Array(256);
let parentCost = Number.POSITIVE_INFINITY;
let parentNet;
let mut;

function setup() {
  createCanvas(400, 400);
  parentNet = new SWNet4(256, 4);
  mut = new Mutator(5, 25, parentNet.scale);
  c1 = color("gold");
  for (let i = 0; i < 8; i++) {
    ex[i] = new Float32Array(256);
  }
  for (let i = 0; i < 127; i++) {
    // Training data
    let t = (i * 2 * PI) / 127;
    ex[0][2 * i] = sin(t);
    ex[0][2 * i + 1] = sin(2 * t);
    ex[1][2 * i] = sin(2 * t);
    ex[1][2 * i + 1] = sin(t);
    ex[2][2 * i] = sin(2 * t);
    ex[2][2 * i + 1] = sin(3 * t);
    ex[3][2 * i] = sin(3 * t);
    ex[3][2 * i + 1] = sin(2 * t);
    ex[4][2 * i] = sin(3 * t);
    ex[4][2 * i + 1] = sin(4 * t);
    ex[5][2 * i] = sin(4 * t);
    ex[5][2 * i + 1] = sin(3 * t);
    ex[6][2 * i] = sin(2 * t);
    ex[6][2 * i + 1] = sin(5 * t);
    ex[7][2 * i] = sin(5 * t);
    ex[7][2 * i + 1] = sin(2 * t);
  }
  textSize(16);
}

function draw() {
  background(0);
  loadPixels();
  for (let i = 0; i < 100; i++) {
    mut.mutate(parentNet.params);
    let cost = 0;
    for (let j = 0; j < 8; j++) {
      parentNet.recall(work, ex[j]);
      cost += costL2(work, ex[j]); // autoassociate
    }
    if (cost < parentCost) {
      parentCost = cost;
    } else {
      mut.undo(parentNet.params);
    }
  }
  fill(c1);
  for (let i = 0; i < 8; i++) {
    for (let j = 0; j < 255; j += 2) {
      set(25 + i * 40 + 18 * ex[i][j], 44 + 18 * ex[i][j + 1]);
    }
  }
  for (let i = 0; i < 8; i++) {
    parentNet.recall(work, ex[i]);
    for (let j = 0; j < 255; j += 2) {
      set(25 + i * 40 + 18 * work[j], 104 + 18 * work[j + 1]);
    }
  }
  updatePixels();
  text("Training Data", 5, 20);
  text("Autoassociative recall", 5, 80);
  text("Iterations: " + frameCount, 5, 150);
  text("Cost: " + parentCost.toFixed(3), 5, 170);
}
